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ABSTRACT 

Daddi et al. have recently reported strong clustering of a population of red galaxies 
at z~3 in the Hubble Deep Field-South. Fitting the observed angular clustering with 
a power law of index —0.8, they infer a comoving correlation length ro ~ 8/i _1 Mpc; 
for a standard cosmology, this ro would imply that the red galaxies reside in rare, 
M > 10 h M@ halos, with each halo hosting ~ 100 galaxies to match the number 
density of the population. Using the framework of the halo occupation distribution 
(HOD) in a ACDM universe, we show that the Daddi et al. data can be adequately 
reproduced by less surprising models, e.g., models with galaxies residing in halos of mass 
M > M min = 6.3 x lO n /i _1 M and a mean occupation iV av g(M) = 1.4(M/M min ) a45 
above this cutoff. The resultant correlation functions do not follow a strict power 
law, showing instead a clear transition from the one-halo-dominated regime, where 
the two galaxies of each pair reside in the same dark matter halo, to the two-halo- 
dominated regime, where the two galaxies of each pair are from different halos. The 
observed high-amplitude data points lie in the one-halo-dominated regime, so these 
HOD models are able to explain the observations despite having smaller correlation 
lengths, ro ~ 5/i _1 Mpc. HOD parameters are only loosely constrained by the current 
data because of large sample variance and the lack of clustering information on scales 
that probe the two-halo regime. If our explanation of the data is correct, then future 
observations covering a larger area should show that the large scale correlations lie 
below a # -L8 extrapolation of the small scale points. Our models of the current data 
suggest that the red galaxies are somewhat more strongly clustered than UV-selected 
Lyman-break galaxies and have a greater tendency to reside in small groups. 

Subject headings: galaxies: halos - galaxies: high-redshift - large-scale structure of 
universe 



1. Introduction 



Clustering of high-redshift galaxies can provide information about their relation to the under- 
lying mass distribution and their formation mechanisms. Efforts have been made to detect high-z 
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galaxies and to estimate their clustering properties. For example, surveys of z ~ 3 Lyman break 
galaxies (LBGs) (Adelberger et al. 1998, 2003; Steidel et al. 1998) show that these galaxies are 
strongly clustered, with a correlation length ro ~ 4/i~ 1 Mpc. This strong clustering appears to 
be naturally explained by theoretical models, which predict high bias of luminous high-z galaxies 
(Governato et al. 1998; Kauffmann et al. 1999; Cen & Ostriker 2000; Benson et al. 2001; Pearce et 
al. 2001; Somerville et al. 2001; Yoshikawa et al. 2001; Weinberg et al. 2004). Recently, using VLT 
observations, Daddi et al. (2003) have analyzed the clustering properties of K < 24 galaxies in the 
Hubble Deep Field South (HDF-S). They find that a population of red galaxies with J — K > 1.7 in 
the photometric redshift range 2 < z p hot < 4 exhibit remarkably strong clustering, rg ~ 8/i _1 Mpc. 
This paper attempts to interpret these measurements in the framework of halo occupation distri- 
bution (HOD) models (see, e.g., Seljak 2000; Scoccimarro et al. 2001; Berlind k Weinberg 2002 
and references therein). 

Fitting the measured two-point angular correlation function by a power law with an index 
—0.8, which corresponds to a power law real-space two-point correlation function with an index 
— 1.8, Daddi et al. (2003) derive a correlation length as large as ro ~ 8/i _1 Mpc for the red galaxies. 
This strong clustering seems hard to reconcile with conventional models of galaxy bias. For a 
reasonable cosmological model, such as that assumed in the GIF simulation of Jenkins et al. (1998), 
a correlation length of ~ 8/i _1 Mpc corresponds to a linear galaxy bias factor of ~ 5 at z ~ 3. In 
the halo bias model (e.g. Mo & White 1996), this bias factor implies that these red galaxies 
would be hosted by M > 1O 13 M halos. The comoving number density of M > 1O 13 M halos 
is ~ 3 x 10~ 5 /i 3 Mpc~ 3 . To match the comoving number density, ~ 7 x 10 _3 /i 3 Mpc -3 , of the red 
galaxies, there should be more than 200 such galaxies in each halo. Even if we take into account 
the fact that galaxy bias is an average of halo bias weighted by occupation numbers and lower the 
halo mass to M > 5 x 10 12 M Q , the occupation number is still as large as about 70. Based on the 
data, Daddi et al. (2003) speculate that the problem may be caused by the effect of a small scale 
excess in the correlation function. Detailed modeling is necessary to resolve this puzzle. 

For modeling observed galaxy clustering statistics, the framework of the HOD is a powerful 
tool. The HOD describes the relation between the distribution of galaxies and that of the matter 
at the level of individual dark matter halos. It characterizes the probability distribution P(N\M) 
that a halo of mass M contains N galaxies of a given type and specifies the relative spatial and 
velocity distributions of galaxies and dark matter within halos. With an assumed cosmological 
model that determines the halo population, the HOD can be inferred empirically from observed 
galaxy clustering (Peacock & Smith 2000; Marinoni & Hudson 2002; Berlind & Weinberg 2002). 
HOD modeling has been applied to galaxy clustering data from the Two-Degree Field Galaxy 
Redshift Survey (2dFGRS) and the Sloan Digital Sky Survey (SDSS) (see, e.g., van den Bosch et 
al. 2003; Magliocchetti & Porciani 2003; Zehavi et al. 2004). HOD modeling has also been used to 
model the clustering of high-z galaxies, such as LBGs (Bullock, Wechsler, &: Somerville 2002) and 
extremely red objects (EROs) (Moustakas & Somerville 2002). 

In this paper, we will apply HOD modeling to the population of red galaxies at z ~ 3 in Daddi 
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et al. (2003) and try to understand the apparent strong clustering of these galaxies. We describe 
the HOD parameterization and how we analytically calculate the galaxy correlation function in § 2. 
In § 3, we explain what we learn from model fitting to the observational data. We summarize the 
results and give a brief discussion in § 4. 

2. HOD Parameterization and Analytic Calculation of Correlation Function 

Motivated by measurements of the cosmic microwave background (e.g., Netterfield et al. 2002; 
Pryke et al. 2002; Spergel et al. 2003), the abundance of galaxy clusters (e.g., Eke, Cole, & Frenk 
1996), and high redshift supernova observations (e.g., Riess et al. 1998; Perlmutter et al. 1999), 
we adopt a spatially flat ACDM cosmological model with matter density parameter Q m = 0.3 
throughout this paper. For the matter fluctuation power spectrum, we adopt the parameterization 
of Efstathiou, Bond, & White (1992) and assume that the spectral index of the inflationary power 
spectrum is n s = 1, the rms fluctuation (linearly evolved to z = 0) at a scale of 8/i _1 Mpc is as = 0.9, 
and the shape parameter is T = 0.21. The Hubble constant is assumed to be WOh kms^ 1 Mpc -1 
with h = 0.7. 

2.1. HOD parameterization 

To do an analytical calculation of the correlation function of galaxies, we need to parameterize 
the halo occupation distribution. For the functional form of halo occupation number, we adopt a 
simple model similar to that used by Zehavi et al. (2004), which is loosely motivated by results from 
smoothed particle hydrodynamics (SPH) simulations and semi-analytic calculations (see Berlind et 
al. 2003 and references therein). In this model, in halos of mass M > M±, the mean occupation 
number follows a power law, iV av g(M) = (M/M\) a , and in halos of mass M m i n < M < Mi there is 
only a single galaxy that is above the luminosity threshold, i.e., N avg (M) = 1. Given a and Mi, 
M m i n is then fully determined by the number density of galaxies. Since the correlation function is a 
statistic of galaxy pairs, we also need to know the second moment of the occupation number. SPH 
simulations and semi-analytic models predict that the distribution of halo occupation numbers at 
fixed halo mass is much narrower than a Poisson distribution when the occupation is low. Here we 
adopt the so-called nearest-integer distribution for P(iV|iV av g), which states that the occupation 
number for a halo of mass M is one of the two integers bracketing N avg (M), with the relative 
probability determined by having the right mean. Besides this basic model, we will also consider 
some alternatives as discussed in § 3. More detailed discussions of the parameterization of HOD 
models appear in Berlind et al. (2003) and Zehavi et al. (2004). Our parameterization here is quite 
restrictive, but the data are not sufficient to constrain a model with more freedom. 
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2.2. Real Space Correlation Function 

The two-point correlation function of galaxies £ (r) reflects the excess probability over a random 
distribution of finding galaxy pairs with a separation r. From the point view of the halo model, the 
two galaxies of each pair can come from either a single halo or two different halos. Consequently, 
the two-point correlation function £(r) can be decomposed into two components, 

£(r) = [l + £ih(r)]+&h(r), (1) 

where the one-halo term £ih(0 and the two-halo term £211 ( r ) represent contributions by pairs from 
single halos and different halos, respectively. The above expression comes from the fact that the 
total number of galaxy pairs (oc 1 + £(r)) is simply the sum of the number of pairs from single halos 
(oc 1 + £ih(j)) an d that from different halos (oc 1 + £;2h( r ))- The one-halo term and two-halo term 
dominate respectively at small and large separations. 

The one-halo term £m(r-) can be exactly computed in real space through (Berlind & Weinberg 
2002) 

l+* ( r ) — 1 r , Af dn (N(N-l)) M 1 / r \ 

1 + " 2^| J dM dM 2 2R~(M) F [2R- ) ' (2) 

where n g is the mean number density of galaxies, dn/dM is the halo mass function (Sheth & 
Tormen 1999; Jenkins et al. 2001), (N(N — 1))m/2 is the average number of pairs in a halo of mass 
M, and F(r/2R V - U ) is the cumulative radial distribution of galaxy pairs, i.e. the average fraction of 
galaxy pairs in a halo of mass M (virial radius i? v ir) that have separation less than r. The function 
F'(x) depends on the profile of the galaxy distribution p g (r) within the halo. In this paper, we 
assume that there is always a galaxy located at the center of the halo, and others are regarded as 
satellite galaxies. With this assumption of central galaxies, F'(x) is then the pair-number weighted 
average of the central-satellite pair distribution F^(x) and the satellite-satellite pair distribution 
F s ' s (x) (see, e.g., Berlind & Weinberg 2002; Yang, Mo, & van den Bosch 2002), 

M^V(*) = (N - + (3) 

The central-satellite galaxy pair distribution F^(x) is just the normalized radial distribution of 
galaxies (i.e., oc p g (r)r 2 ), and the satellite-satellite galaxy pair distribution F^ s (x) can be derived 
through the convolution of the galaxy distribution profile with itself (see Sheth et al. 2001). We 
will first assume that the galaxy distribution is the same as the dark matter distribution within the 
halo, which follows an NFW profile (Navarro, Frenk, & White 1995, 1996, 1997) truncated at the 
virial radius. The concentration of an NFW profile depends on the halo mass, for which we use the 
relation given by Bullock et al. (2001) after modifying it to be consistent with our halo definition - 
a gravitationally bound structure with overdensity ~ 200. Later in this paper, we will also consider 
a more concentrated galaxy distribution profile. 



The two-halo term is basically the average halo correlation function weighted by the average 
occupation number of galaxies of each halo. The halo correlation function is related to the mass 
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correlation function by the halo bias factor (Mo & White 1996; Jing 1998; Sheth, Mo, & Tormen 
2001). It is convenient to calculate the two-halo term in Fourier space (Seljak 2000; Scoccimarro 
et al. 2001). The two-halo term contribution to the galaxy power spectrum reads 

2 

(4) 

where P m (k) is the mass power spectrum, N avg (M) is the mean occupation number in halos of 
mass M, bh(M) is the halo bias factor, y g (k,M) is the (normalized) Fourier transform of the 
galaxy distribution profile within a halo of mass M, and M max is the upper limit for the integral 
(see below). In the calculation, we adopt the three improvements mentioned in Zehavi et al. (2004). 
First, for P m (k), instead of the linear spectrum as used in previous studies, we use the non-linear 
power spectrum as given by Smith et al. (2003) to account for the non-linear evolution of the 
mass (also see Magliocchetti & Porciani 2003). Second, the halo exclusion effect is approximately 
considered by choosing an appropriate M max : for the two-halo term at separation r, M max is taken 
to be the mass of the halo with virial radius r/2. Third, the scale-dependence of the halo bias 
factor on non- linear scales is incorporated by using an empirical formula from simulations. The 
two-halo term of the correlation function is the Fourier transform of the power spectrum, 

^) = ^[ P ^k)k 2 ^dk. (5) 

The correlation function analytically calculated using the above method agrees fairly well with 
that measured from a mock galaxy catalog generated by populating galaxies according to the same 
HOD into halos identified in iV-body simulations (see Zehavi et al. 2004 and Figure 2 below). 



P m (k) — P m {k) 
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2.3. Angular Correlation Function 

The angular distribution of galaxies is a projection of the three-dimensional distribution. The 
angular correlation function w{6) of galaxies is related to the real-space correlation function through 
Limber's equation (Peebles 1980). In a flat universe, as adopted in this paper, for the small-angle 
limit, Limber's equation has the form 

™ -CT n 2 g {x)S\x) X ±dx J A l e(Vy 2 + ^ 2 , z)dy 
w{6) = — = -J , (6) 

n g (x)S(x)x*dx 



where Ar = r max — r m j n is the radial range of the survey, n g (r) is the average number density of 
galaxies at distance r, and S(r) is the selection function of the sample (all distances are in comoving 
units). The sample selection function S(r) can be derived from the observed redshift distribution 
if the sample is large enough, but the 49 galaxy sample of Daddi et al. (2003) is not large enough 
to allow a precise measurement. However, it turns out that the basic result of this paper is not 
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sensitive to the form of the selection function. We therefore simply assume that S(r) is constant 
over the redshift range z = 2 to z = 4, which defines r m i n and r a 



max- 



In practice, the angular correlation function is estimated by comparing the observed pair 
numbers in an angular separation bin with those from a random sample of the same geometry. 
The widely used estimator proposed by Landy & Szalay (1993) estimates the angular correlation 
function as 

DD - 2DR + RR 

MO) = > (7) 

where DD, DR, and RR represent number counts of data-data (galaxy-galaxy) pairs, data-random 
(galaxy-random) pairs, and random-random pairs, respectively, in the angular bin around 9. Each 
of these number counts are normalized so that the summation over all 9 is unity (i.e., the number 
of pairs in each angular bin is divided by the total number of pairs in the field). The estimated 
angular correlation function Wb(9) is subject to a statistical bias that leads to systematically lower 
values than the true angular correlation function w(9), Wb(9) = w(9) — I.C., where 

I.C. = ^2 J J u>(0)dftidft 2 (8) 

is the integral constraint (Groth & Peebles 1977). Since Wb(9) is the quantity directly measured 
from the observation, it is more appropriate to try to fit Wb{9) than to fit w(9). To convert the 
analytically predicted w(9) to u>b(9), we use the random-random sample to calculate the integral 
constraint expected for the model w(9) (see, e.g., Roche et al. 1999), 

where N TT (9) is the count of random-random pairs in the angular bin around 9. We only need to 
estimate f{9) = N„{9) /T,N„{9) once from a random sample that has the same geometry as the 
observation. We generate 100 such random samples with 5,000 points in each and take the average 
f(9) as the estimate. 



3. Fitting the Observations 

The angular clustering data we are interested in are for a population of selected galaxies 
{K < 24) at 2 < Zphot < 4 with J — K > 1.7. Details about this sample can be found in Daddi et al. 
(2003). The sample includes 49 galaxies found within a field of view of ~4 arcmin 2 . The comoving 
number density of these galaxies is ~ 7.1 x 10~ 3 /i 3 Mpc -3 . Assuming the angular correlation function 
to be a power law with an index —0.8, Daddi et al. (2003) find its amplitude at 1° to be 39.1 ± 10.2, 
which corresponds to a real space correlation length of 8.3 ± 1.2/i _1 Mpc (comoving). 

We now fit the data (kindly provided in electronic form by E. Daddi) using the model in § 2. For 
a given assumption about P(N\N avg ), e.g., a nearest-integer or Poisson distribution, the analytic 
angular correlation model discussed in § 2.3 has two free parameters: Mi, which determines the 
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amplitude of N avg (M), and a, which is the slope of iV av g(M) at high halo masses. M m i n is fixed by 
the mean number density of galaxies in the sample. We thus perform a two-parameter \ 2 fit to the 
data. The observational error bars reported by Daddi et al. (2003) are used in the calculation of 
X 2 and are assumed to be uncorrelated. With these data and error bars, we find that the two free 
HOD parameters are highly degenerate and that they can be only loosely constrained individually. 
For example, with the nearest-integer distribution, Mi is in the range 4-30x lO lo /i _1 M and a is 
in the range 0.2-0.5. If we assume a Poisson distribution for P(N\N lW g), then the resultant a is 
unrealistically large (~ 3), with large uncertainty. With the Poisson distribution, if we change the 
functional form of N avg (M) to be a power law with a low mass cutoff, a can be in a reasonable 
range but still with large uncertainty. Through studying subhalos in high-resolution dissipationless 
simulations, Kravtsov et al. (2004) propose an HOD form, which separates contributions from 
central and satellite galaxies. The mean occupation function of central galaxies is a step function, 
while the distribution of satellite galaxies can be approximated by a Poisson distribution with the 
mean following a power law. The resultant shape of the mean occupation function and scatter 
around the mean are somewhat similar to our basic model. We also try this HOD form and find 
results and uncertainties similar to those of our basic model. 

We illustrate the looseness of the HOD parameter constraints in Figure 1 by showing the 
results of different parameter combinations that lead to similar real-space correlations and angular 
correlations for the nearest-integer case. Note that since M m - m > M\ is derived from the fit, the 
resultant A r avg (M) is equivalent to a case in which iV aV g(M) is a power law with a low-mass cutoff, 
with no "flat" portion at M < M\. The result of M m [ n > M\ mimics the case of local giant elliptical 
galaxies and z ~ 1 EROs as modeled by Moustakas & Somerville (2002), a point discussed further 
at the end of this section. There is a strong break in the modeled correlation function between 
the one-halo and two-halo regimes. The sharpness of this break is somewhat exaggerated by the 
approximate nature of our correction for halo exclusion. However, the angular correlation function 
is less affected by this approximation because it is a projection of the real-space correlation, and 
we show later that the approximate treatment of halo exclusion has a negligible effect on our Wb(9) 
modeling here. 

The reduced x 2 (7 degrees of freedom, 9 data points, and 2 free parameters) from either the 
nearest-integer model or the power law Poisson model is about 1.8, which does not seem to be 
a good model fit. The field of view of the survey is less than 4 arcmin 2 , and the total number 
of galaxies is only about 50. Thus, as noted by Daddi et al. (2003), the sample variance may be 
large, and error bars based only on the finite number of objects (as used above) may be too small. 
We therefore attempt to make more realistic error estimates by populating halos from the GIF 
simulation (Jenkins et al. 1998) with galaxies. We use the halo population from the GIF simulation 
output at z = 2.97 and proceed as follows. First, each halo is assumed to have a truncated NFW 
density profile with the same concentration-mass relation used in the analytic model. We then 
populate galaxies according to the halo occupation distribution from the \ 2 fitting and generate 
a mock galaxy catalog. Next, we randomly extract 10 slices along one direction from the mock 
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catalog, with the cross-section of each slice having the same size and geometry as the observation. 
These 10 slices are checked at selection to make sure that they do not overlap (even partially) with 
each other. The 10 slices are stacked together to approximate the radial extent of the survey in 
comoving distance. The projection of the stacked slice thus represents one "observation." Finally, 
we estimate Wb(9) for this observation in the same angular bins as the real data using the technique 
of Landy & Szalay (1993) (Equation 7). The data-random and random-random terms are averaged 
over 100 random realizations, and each random sample realization has 5000 points. Altogether we 
make 100 observations and estimate Wb(0) for each one. 

The result is shown in the top left panel of Figure 2. The central solid line is the average 
over the 100 observations, which agrees with the model prediction (dot-dashed line) as expected 
(and verifying that our analytic approximation is accurate enough for our purposes). The dashed 
lines above and below the solid line represent the la scatter of the 100 observations. The estimated 
angular correlation Wb(0) for an individual observation is very uncertain and may even not decrease 
monotonically with 6 (as is the case for the real data points). Compared with the scatter derived 
here, the observational error bars are apparently underestimated by a factor of about 1.5. If we 
take the mock catalog scatter as true error bars, then our model fit is acceptable. However, the 
uncertainties in HOD parameters were large even with the original error bars, so we do not wish 
to place much weight on the particular values that emerge in the best fit. Rather, we wish to 
use our HOD models as a general guide in understanding the implications of the data. Perhaps 
the most important lesson is that the observed angular correlation signals are dominated by the 
one-halo term, where the two galaxies of each pair are from one single halo. This can be seen in 
the top left panel of Figure 2, where the dotted line shows the two-halo term. The contribution 
from the two halo term becomes comparable to that from the one-halo term only on angular scales 
greater than ~ 0°005 (also see the bottom left panel of Figure 1). As mentioned in Daddi et al. 
(2003), the estimated angular correlation at the smallest angular scale is mainly from a few triplets. 
The redshift distribution of the galaxy sample provides further evidence. There are many spiky 
structures in the redshift distribution (Daddi et al. 2003). Since the largest projected separation in 
the field of view is about 3/i _1 Mpc, galaxies within the same spike are most likely to be physically 
close, and thus they have a high probability of being located in the same halo. 

Domination of the signal by the one-halo term has several implications. The HOD model 
generically predicts that the correlation function is not strictly a power law (see Berlind & Weinberg 
2002). Instead, there should be a transition region from the one-halo-dominated regime at small 
scales to the two-halo-dominated regime at large scales. For the real-space correlation function, 
such a transition happens around the virial radius of the largest halos, 2 — 3/i _1 Mpc at z = 0. 
Recent results from SDSS have revealed a statistically significant departure from a power law in 
the two-point correlation function, which can be well explained within the framework of the HOD 
(Zehavi et al. 2004). Two-point correlation functions measured from other surveys, such as 2dFGRS 
and APM, also show such a departure (e.g. Hawkins et al. 2003; Padilla & Baugh 2003). 

The two-point correlation function predicted by the model that fits the angular correlation 
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function in this paper shows a prominent departure from a power law (Figures 1 and 2), which 
is also reflected in the predicted angular correlation function (Figure 1). [Note that the excellent 
agreement between the numerical and analytic calculations of Wb(0) demonstrates that any artifacts 
of our approximate treatment of halo exclusion are negligible in comparison to the observational 
error bars.] Since the one-halo term is related to the distribution of galaxies within halos and the 
two-halo term is mostly determined by the halo-halo distribution, the amplitude and slope of the 
two terms may differ from each other substantially. One should therefore be cautious about inferring 
the correlation length by assuming a pure power law in the correlation function. In Daddi et al. 
(2003), a power law with an index —1.8 for the real-space correlation function, corresponding to 
an index of —0.8 for the angular correlation function, is assumed, and a high correlation amplitude 
(correlation length ~ 8/i _1 Mpc) is found. This strong correlation is unlikely to be related to the real 
correlation between halos, since, as we show here, the statistically significant signal is dominated 
by galaxy pairs within halos. The distribution of galaxies within halos does not tell how galaxies 
cluster on large scales, and the correlation length is overestimated because of the high amplitude 
of the one-halo term. In fact, from the fitting model, the correlation length where £(r) = 1 can 
be as low as ~ 5/i _1 Mpc (Figure 2), which is in a good agreement with the result of Kravtsov 
et al. (2004) for a subhalo sample of comparable number density. The mystery about the strong 
clustering in this particular sample then disappears. Our explanation of this mystery is, in some 
sense, a quantitative version of the speculation of Daddi et al. (2003) that the strong clustering 
signal is an effect of "excess" small-scale clustering. 

We note that although the fit to the data can be regarded as acceptable, the third and the 
fourth data points are well below the prediction. This may be of no significance considering the 
large sample variance. Nevertheless, it is interesting to ask what the cause may be if this discrep- 
ancy is real. The low amplitude of these two data points may be a hint that the one-halo term 
drops faster than in our model, which means that the distribution of galaxies within halos is more 
concentrated than the NFW profile we use. As an alternative model, we first doubled the con- 
centration parameter of the galaxy distribution profile within each halo, thus making the galaxies 
more centrally concentrated than the dark matter, but this change is not adequate to match the 
observed w\,(Q) in the third and fourth angular bins. It thus implies that the distribution profile 
of galaxies is steeper than the NFW profile. As a more extreme alternative, we considered an r~ 3 
profile, with a flat core at r < 0.1i? v i r to make the pair distribution finite. The bottom panels of 
Figure 2 show the resultant model fitting and the sample variance estimated from mock catalogs 
generated using halos in the GIF simulation. The steeper galaxy distribution yields a better fit to 
the third and fourth data points. As before, HOD parameters remain poorly constrained. With 
the current sample size, the preference for r -3 profiles over NFW profiles is not highly significant, 
but the low amplitude of the third and fourth data points could be a hint that the red galaxies in 
the sample of Daddi et al. (2003) are centrally concentrated within their parent halos, analogous 
to the morphological segregation observed in present-day clusters (e.g., Oemler 1974; Melnick & 
Sargent 1977; Dressier 1980; Adami, Biviano, & Mazure 1998). 
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Although HOD parameters are loosely constrained, the cutoff mass M m j n in all the models 
shown in the figures roughly remains constant, ~ 6 x 10 11 /i _1 Mq. The approximate constancy of 
M m i n mainly comes from the constraint of the galaxy number density and the steep drop of the halo 
number density toward higher halo masses. For example, the cumulative number density of halos 
drops from ~ 5 x 10~ 3 /r 3 Mpc~ 3 to ~ 3 x 10 _3 /i 3 Mpc~ 3 as the minimum halo mass changes from 
5 x 10 11 /i _1 M Q to 7x 1O 1:l /i _1 M0. With the galaxy number density fixed, a large change in M m - m is 
hard to compensate with changes in other HOD parameters. Although the sharp cutoff in N avg (M) 
that we have assumed in this paper is an idealization, the derived value of M m i n should still give an 
approximate indication of the characteristic minimum masses of halos that host the red galaxies. 
In our successful models, the mean occupation number at M m j n is above 1 (i.e., M m \ n > Mi). 
This suggests that the red galaxies arise preferentially in groups and clusters (see Moustakas & 
Somerville 2002), which may be a signature of an environmental effect on color. However, since 
M m i n < M2, where M2 is the mass of the halo that on average contains two red galaxies, there are 
still single-occupancy halos as the nearest-integer distribution is taken into account. For example, 
in the model with a = 0.45 and M min = 6.3 x lO n /i _1 M , M 2 = 1.4 x 1O 12 /i _1 M , and about 11% 
of the galaxies are the sole occupants of their halos. 

There are several hints that the red galaxies of the Daddi et al. (2003) sample have clus- 
tering properties different from those of the UV-selected galaxies (e.g., LBGs) at the same red- 
shift. The first hint comes from the M m \ n > M\ result mentioned above. Using similar kinds of 
HOD models (although assuming a pure power law with a low mass cutoff, with no single occu- 
pancy "plateau"), Bullock, Wechsler, & Somerville (2002) and Moustakas k Somerville (2002) find 
M min ~ lO lo /i _1 M < Mi for LBGs at z ~ 3, implying that most LBGs are the sole occupants of 
their parent halos. By contrast, our model fits imply that most red galaxies reside in groups of two 
or more, similar to the results of Moustakas & Somerville (2002) for local giant elliptical galaxies 
and z ~ 1 EROs (for which they find a trend of M m i n > Mi). A second hint is from the correlation 
length itself, which is still ~ 5/i _1 Mpc in our models. The correlation length of UV-selected LBGs 
in the spectroscopic sample of Adelberger et al. (2003) is only about 4/i _1 Mpc, and it appears to 
decrease for samples of lower luminosity threshold and higher space density (Giavalisco & Dickinson 
2001). A final hint comes from the behavior of clustering on small scales. Porciani & Giavalisco 
(2002) find that the angular correlation function of z ~ 3 LBGs drops at separations of less than 
30", while the angular clustering of the red galaxies seems, if anything, to be exceptionally strong 
at the smallest angular scales. To illustrate these differences between the red galaxies and the 
UV-selected LBGs, we show in Figure 1 what the HOD and correlation functions would look like 
for an LBG sample that has the same number density as the red galaxies. For this purpose, we 
start from the HOD model result of z ~ 3 LBGs by Moustakas & Somerville (2002), which is 
consistent with that of Bullock, Wechsler, & Somerville (2002), and change their HOD parameters 
a little bit to match the number density here. The mean occupation function has a power-law 
form iVavg(M) = (M/Mi) a with a low-mass cutoff M min . We adopt M min = 1.4 x lO lo /i. _1 M , 
Mi = 4.0 x 10 12 /i~ 1 M Q , a = 0.8, and a nearest-integer distribution. The dashed line in the top 
right panel of Figure 1 shows this mean occupation function, where we can clearly see that unlike 
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the red galaxies, relatively fewer LBGs reside in groups. This leads to a lower small-scale clustering 
amplitude and a lower correlation length with respect to the red galaxies (see the dashed curves 
in the left panels of Figure 1). For a more consistent comparison between the red galaxies and the 
LBGs, we need detailed HOD modeling of the LBGs, which is out of the scope of this paper. The 
exercise here is to simply illustrate the differences in the clustering properties of the red galaxies 
and the LBGs, as noticed by Roche et al. (2002), Roche, Dunlop, & Almaini (2003), and Daddi et 
al. (2003) from angular clustering measurements. 

4. Conclusion and Discussion 

In this paper, we have presented an HOD model of the observed strong clustering of a pop- 
ulation of red galaxies at z ~ 3 analyzed by Daddi et al. (2003). Fitting the angular correlation 
data by assuming the real-space correlation to be a power law with the form (r/ro)" 1 ' 8 , Daddi et 
al. (2003) find the correlation length ro ~ 8/i -1 Mpc, which would imply that galaxies reside in rare 
halos with M > 1O 13 /i _1 M and which would require a very large occupation number in each halo 
to account for the observed galaxy number density. Our HOD modeling shows that the angular 
clustering data can be explained by a less surprising model, e.g., with a cutoff at 6.3 x 10 n /t M 
and mean galaxy occupation number iV avg (M) = 1.4(M/6.3 x W u h~ 1 M Q ) 0A5 above this cutoff. 
Artificial galaxy catalogs constructed with this HOD show that sample variance increases error 
bars by ~ 50% over those estimated by Daddi et al. (2003), which (as they noted) did not take 
sample variance into account. 

There is degeneracy between HOD parameters M\ and a. However, the characteristic minimum 
mass of halos that can host the red galaxies seems to be around 6 x 10 11 /i~ 1 M Q . Results from our 
modeling suggest that the red galaxies are a different population from LBGs. 

HOD parameters are not tightly constrained by the data, but in all cases the significantly 
non-zero points from Daddi et al. (2003) are in a regime dominated by pairs within single halos. 
The amplitude of the correlation function in the two-halo regime is below an r~ L8 power law 
extrapolation of that of the one-halo regime, which is why lower mass halos are acceptable. Thus, 
if our explanation is correct, surveys with larger area should show weaker correlations than this 
r extrapolation. The correlation length predicted by our model can be as low as ~ 
a prediction that can be tested by larger area surveys. 

Obtaining good constraints on HOD parameters will require samples large enough to accu- 
rately probe the two-halo regime as well as the one-halo regime. Spectroscopic surveys are also of 
importance since they allow one to measure galaxy clustering in redshift space (in addition to more 
accurate measurements of the projected clustering). With a good understanding of the velocity field 
of halos, the clustering in redshift space would at least provide a consistency check for the HOD 
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model. With wider angle space-based surveys such as GOODS 1 and ambitious infrared follow-up 
programs like the FIRES project (Franx et al. 2000; Daddi et al. 2003), the necessary data should 
become available in the next several years. This will provide valuable constraints on the host halos 
of red high-z galaxies and clues to their formation histories. 

We thank Emanuele Daddi not only for providing us the observation geometry and the angular 
clustering data used in this paper but also for his useful comments and suggestions. We also thank 
David Weinberg for helpful discussions, for his encouragement on this work, and for his valuable 
comments that improved the paper. This work was supported by NSF grant AST 00-98584 and by 
a Presidential Fellowship from the Graduate School of the Ohio State University. 
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Fig. 1. — Illustration of the loose constraints on individual HOD parameters for the nearest-integer 
model. Three cases of parameter combinations are shown: (M m i n , M±,a) = (6.5 x 10 n ,4.5 x 
10 10 ,0.22), (6.2 x 10 n ,2.0 x 10 n ,0.36), and (6.3 x 10 n ,2.9 x 10 n ,0.45), where masses are in units 
of h~ 1 M & . The top right panel shows the corresponding mean occupation number as a function 
of halo mass for the three cases. The top left panel, the bottom left panel, and the bottom right 
panels are for the real space two-point correlation function, the angular correlation function, and the 
measured angular correlation function (i.e., with the integral constraint subtracted), respectively. 
The dotted lines show contributions from one-halo pairs and two-halo pairs. Data points with error 
bars in the bottom right panel are from Daddi et al. (2003). The dashed line in the top right panel 
illustrates the mean occupation function of the LBGs, and the dashed curves in the left panels are 
corresponding correlation functions (see the text). 
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Fig. 2. — Fitting results, sample variance from mock catalogs, and comparison for different 
assumptions about the galaxy distribution profile within halos. The nearest-integer distribu- 
tion is used. Top : Galaxies distributed according to the NFW profile and (M m [ n , M±,a) = 
(6.3 x lO n /i _1 M ,2.9 x lO n /i _1 M ,O.45). Bottom : Galaxies follow an r~ 3 distribution pro- 
file and (M min ,Mi,a) = (5.8 x 10 11 /i" 1 M Q , 2.7 x 1O 11 /i" 1 M , 0.38). Left panels show the measured 
angular correlations Wb(9). The solid, dashed, and dotted lines are the mean, la scatter about the 
mean, and the two-halo pair contribution, respectively, measured from mock galaxy catalogs gen- 
erated through populating z = 2.97 halos from the GIF simulation (see the text). The dot-dashed 
line is the analytic prediction of Wb(0). Data points with error bars are from Daddi et al. (2003). 
The right panels show the corresponding real-space two-point correlation functions. The dotted 
lines are the one-halo and two-halo terms. Arrows indicate ro where £(ro) = l. Two power law 
curves, (r/ro)~ L8 (dashed curve) and (r/8.3/i _1 Mpc) -1 ' 8 (dot-dashed curve), are also plotted for 
comparison. 



